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Clusters of the two-dimensionnal t-J model with 2 holes and up to 26 sites are 
diagonalized using a Lanczos algorithm. The behaviour of the binding energy with 
system size suggests the existence of a finite critical value of J above which binding 
occurs in the bulk. Only the d-wave pair field operator acting on the Heisenberg GS 
has a finite overlap with the 2 hole ground state for all the clusters considered. The 
related spectral function associated with the propagation of a d-wave (singlet) pair 
of holes in the antiferromagnetic background is calculated. The quasiparticle peak 
at the bottom of the spectrum as well as some structure appearing above the peak 
survive with increasing cluster size. Although no simple scaling law was found for the 
quasiparticle weight the data strongly suggest that this weight is finite in the bulk 
limit and is roughly proportional to the antiferromagnetic coupling J (for J < 1). 
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I. INTRODUCTION 



The possibility that high-temperature superconductivity in the quasi-two dimensional 
cuprates would be a new phenomenon based on purely electronic interactions has moti- 
vated a huge theoretical effort to better understand the physics of two dimensional strongly 
correlated electrons. The Hubbard model and its large coupling version the t-J model are 
one of the simplest theoretical models used to describe the low-energy excitations of the 
copper oxide planes. Although it is difficult to develop a satisfactory perturbative analysis 
in the case of strongly repulsive interactions, numerical methods can often provide useful in- 
formations on both static and dynamical properties of these systems. Exact diagonalization 
(ED) studies do not suffer from random noise problems like stochastic Monte Carlo methods 
in particular in the vicinity of the antiferromagnetic Mott insulator phase. ED are therefore 
well adapted to the calculations of dynamical correlation functions at small doping density. 
However, so far ED have been restricted to fairly small clusters (typically 4x4). Recently, 
an attempt was made to handle larger sizes (up to 26 sites) in order to perform a finite size 
analysis of the data |2|-44| in the close vicinity of the insulator magnetic phase. Although 
this analysis was restricted to the special limit of a single hole in the antiferromagnet (AF) 
it gives some hints about the hole propagation at small but finite doping fraction []2],|3| in the 
region where the antiferromagnetic correlation length is still larger than the spatial extention 
of the hole wavefunction. This work strongly suggested that the quasiparticle (QP) peak 
seen at the bottom of the spectral function survive in the thermodynamic limit || . However 
this analysis neglects the role of a possible hole-hole effective attraction that might occur 
for more than a single hole in the AF. The purpose of the present work is then to study 
the case of 2 holes ie the simplest case that nevertheless include the effect of the effective 
attraction between the holes. The issue of binding is of great importance in the search for 
superconductivity in models of strongly correlated fermions. By a similar finite size analysis 
I shall attempt to give an estimate of the binding energy and of the pair spectral function 
in the thermodynamic limit. Partial results for the binding energy and the QP weight have 
been reported elsewhere ||. The emphasis is put here on the energy dependence of the 
complete spectral function ^^(u;) and on its behavior with system size. 



In standard notations the t-J Hamiltonian reads 

H = JZUSfSfe- irnni+r) + l(S?S7 +g + <?r <?+_.)} 

-tJ2(c\,*c i+la + h.c.). (1.1) 

The first term is the usual antiferromagnetic coupling between spins which splits into a 
diagonal part and a spin flip part. The second term describes the hopping of the holes (ie 
where spins are missing). The sum over i, e is restricted to nearest neighbor bonds along 
x and y on a 2D square lattice, c\ a = (1 — rii- a )ci^ is the hole creation operator and 
rii = J2a c l cr c i,cr- So far t will set the energy scale. 

The Hamiltonian ( |1 . 1| ) with 0, 1 and 2 holes is diagonalized on cluster of increasing sizes 
by a standart Lanczos procedure. It is important in this analysis that all the cluster have 
the same shape namely a square. One should expect in this case a smoother extrapolation 
to the thermodynamic limit. If the (integer) coordinates of the cluster periods are (n,m) 
and (— m, n) the cluster contains N = n 2 + m? sites. The following analysis is based on 
clusters of sizes 4x4, \/l8 x Vl~8, a/20 x and y/26 x y^26. In order to be able to 
handle the largest size (26 sites) use of complete translation and rotation symmetries (for 
and 2 holes) as well as time reversal (for an even number of fermions) is necessary and the 
hamiltonian is diagonalized in each symmetry sector. The size of the Hilbert space increases 
exponantially fast with the cluster size. In the relevant symmetry sector (for discussion about 
GS symmetry see Sec. 2) the Hilbert space for 2 holes contains after symmetry reduction 
3 332, 13 858, 58 274 and 4 229 236 states for N = 16, 18, 20 and 26 sites respectively. 

II. GROUND STATE ENERGY AND BINDING 

Symmetry properties of the Heisenberg GS | \I/q > and of the two-hole (singlet) GS 
| "&o~ 2 > are summarized for various cluster sizes in Tables 1A and IB respectively. The 
symmetries of the half filled GS are exactly known || but they depend on the convention 
to represent the spins. In the following I use the fermion representation including at half 
filling. | ^Iq > and | ^q~ 2 > are both invariant under lattice translations (total momentum 
K = 0) and even (odd) under spin reversal for 4p (4p+2) fermions. This last feature is a 
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simple consequence of the singlet nature of these states. Although, at half filling, in the spin 
language the GS is always invariant under a 90° rotation around any lattice site (s-wave) 
|J, fermion reordering can lead to an extra minus sign, eg for iV = 4 or iV = 20 the half 
filled GS becomes d-wave in the fermion convention (see Table 1A). It is worth noticing that 
in such a case the 2 hole GS is s-wave (see Table IB). However in most cases the half filled 
GS is s-wave and the 2 hole GS is d-wave. 

The physical idea of an effective attraction between the holes can be simply understood 
in the Ising limit of Hamiltonian (|1.1|) ie when the spin flip term is neglected. In this case 
the half filled GS is the simple classical Neel state and obviously two holes added in this 
background can minimize the magnetic energy cost by sitting on nearest neighbor sites. 
Although this crude explaination neglects the delocalisation energy of the holes (t term) this 
picture was indeed shown numerically to be correct provided that J > 0.18 |$). The situation 
is not so clear in the fully quantum case ( |1 . 1|) where spin fluctuations are included. Numerical 
calculations of the binding energy on the 4 x 4 ]9}4T2[, -\/l8 x @ and x |TJ] 



clusters suggested that binding still occur but above a slightly larger critical value of J as 
expected. From a recent calculation of the hole density correlation functions in a restricted 
Hilbert space on clusters with up to 26 sites Prelovsek and Zotos quote a critical value of J 
of order 0.2. However this technique gives a poor accuracy for the GS energy and can not be 
applied to calculate directly the binding energy (see later). On the other hand, Boninsegni 
and Manousakis [[U| using Green Function Monte Carlo simulations carried out on large 
clusters at J > 0.4 and extrapolated to the bulk limit and smaller J/t ratios report a critical 
value of J of order 0.28. 

The GS energies of a single hole Eq -1 — Eq as well as the corresponding GS quantum 
numbers have been reported elsewhere [Q. The two hole GS energies Eq~ 2 — Eq are listed 
in Table 2 for various cluster sizes and J values. I note that the size dependence is weak. 
The data of the 26 site cluster confirm the J u dependence, v ~ 0.9, reported earlier for the 



4x4 cluster Wl 



The binding energy is defined quantitatively by 

A 2h = E»- 2 + E»-2E»- 1 - (2.1) 



Binding occurs when A 2 h < in the bulk. In this case | A 2 h | gives the magnitude of 
the attractive potential between the two holes. The simple (crude) arguments stated above 
would give an attraction of the order of a fraction of J. On the contrary, if binding does not 
take place one would expect A 2 h — > when N — > oo. 

A 2 h vs J is shown in Fig. la and the corresponding data are listed in Table 3. For 
N = 26 Fig. la shows a significant increase of the critical value of J above which binding 
occurs, J c ~ 0.125 for N = 26. In Fig. 2a the same data are plotted vs 1/JV. Although 
A 2 h is a small number corresponding to the difference between two quantities of the same 
magnitude its behavior with system size is rather smooth (although not monotonic). Note 
that the jump of the single hole GS momentum [p|J17 between the different clusters may be 



responsible for a small systematic error. Fig. 2a shows unambiguously that binding occurs 
for, let say, J > 0.5. However the strong size dependence at smaller J does not unable to give 
any accurate estimation of the critical value of J for the bulk. I simply note that even for the 
largest size considered here A 2 h ~ — J while one would expect a somewhat smaller (absolute) 
value in the bulk. This means that these numbers might still be far from the expected ones 
in the thermodynamic limit and it is difficult then to perform any extrapolations. On the 
basis of these data I can nevertheless speculate that the critical value of J lies in the range 
0.3 - 0.5. This is of central importance since binding of holes is a necessary condition (not 
sufficient) for the appearance of superconductivity. 

In this work I have not studied the transition (with increasing J) towards phase separation 



19f1 ie when the holes tend to cluster in a separate region of space and separate from the spins 
rather than bind by pairs (if long range Coulomb repulsion is neglected). In fact, according 
to previous work |T0|JT5[| the hole derealization energy prevents the transition towards phase 
separation to occur at the same critical value of J at which binding first appears but rather 
leads to phase separation at a larger value of J. Therefore, in this scenario, there exists 
a finite range of J where holes do form pairs but do not phase separate from the spins. 
Further work is clearly needed and diagonalization of 4 holes on 26 sites would certainly 
help to clarify this issue. This is left for future work. 



III. QUASIPARTICLE WEIGHT AND SPECTRAL FUNCTION 



Let me now consider dynamical correlations. To study the propagation of a pair of holes 
in the underlying AF background it is useful to define the pair spectral function, 

AM = £ I «- 2 I At I <> f S(u + E» - EZ- 2 ). (3.1) 

m 

As defined above | ^q) is the antiferromagnetic GS at half filling (N spins on N sites). 
The sum is performed over a complete set of eigenstates {\l/^~ 2 } of the two hole sector 
(N-2 spins), with corresponding energies E^~ 2 . In the numerical results below, for plotting 
purposes, the ^-functions are replaced by sharp Lorentzians of small width e ~ 0.02. A* is 
the usual pair creation operator, 



N % 

where the form factor J-{e) (= ±1) can have s- or d-wave symmetry. As seen previously, 
it is remarkable that the half filled and two hole GS have systematically opposite quantum 
numbers (ie ±1) under a 90° rotation (around a lattice site). Therefore it is clear that only 
the d-wave pair operator (ie ^-"(±x) = 1 and ^(iy) = —1) has a non-zero matrix element 
between the two GS. For this reason, I shall restrict myself to this choice so far. Note that 
the pair operator ( |3.2j ) is odd under spin reversal and is invariant under lattice translations 
so that, in addition to its d-wave nature, it has precisely the right translation and spin 
symmetries (it is a singlet) to connect the two GS. In the right hand side of Eq. (|3.1|) the 
sum is then restricted to the states of the same symmetry class as the 2 hole GS. 



The calculation of the spectral function (|3.1| ) is done in three steps. First the Heisenberg 



GS is calculated. In practice the number of Lanczos iterations required to obtain a good 
accuracy of the GS increases slowly with the system size (crudely it is proportional to the 
system size). Typically, for the Heisenberg GS I have reached a high accuracy after 60, 150, 
400 and 600 iterations for N = 16, 18, 20 and 26 respectively. Note that the GS energy 
converges in much less iterations than the corresponding wavefunction, in any case less than 
50 iterations. In the second step the pair operator is applied on the Heisenberg GS. Lastly the 
vector generated in the second step is used as the initial vector of a second Lanczos run in the 



2 hole subspace. The spectral function can be generated by a continued fraction expension 



that can be truncated after N it iterations I have carefully studied the convergence of 



the spectral function as a function of Nu- Fig. 3 shows the behavior of A^h for increasing Na 
in the case of the 26 site cluster. At small Na the spectral weight is distributed on a small 
number of (^-functions. It is important to note that the low energy part of the spectrum is 
the first to converge when Nu is increased. Fig. 3 shows that a complete convergence on 
the whole physically relevant energy range can be obtained after a rather small number of 
iterations much smaller than the size of the Hilbert space. Empirically I have noticed that 
the value of Nu required to obtain a good convergence does not grow much faster than, let 
say, the system size N. For example, N it = 200, 400, 500 and 450 iterations are needed for 
16, 18, 20 and 26 sites. 

The spectral function is plotted in Fig. 4 for increasing cluster size and for J = 0.1, 0.3 
and 1. I observe that most of the features of the 4x4 cluster |2(J survive when the cluster 
size is increased; (i) a QP peak lies at the bottom of the band and (ii) a significant amount 
of spectral weight is spread at higher energy over a range of a few t (for details see discussion 
below). This is very similar to the single hole spectral function I note that the few 

peaks seen in the case of the smallest size (16 sites) are smeared out when N is increased and 
eventually merge into an almost perfect continuous background for N=26 except for J = 1. 
It is also remarkable that the various local maxima in the background (at uj ~ 0.8 in Fig. 4d 
or u ~ 0.5 , 2 and 4 in Fig. 4h) are located at almost the same frequencies for the various 
sizes N considered. However it is not clear whether the secondary peaks (of finite width) 
observed immediately above the QP peak in Figs. 4h and 41 are still finite size effects or 
whether they are related to any physical low energy resonances (see later for discussion). It 
should be noted that, similarly to the single hole case ||, Figs. 4 suggest that no real gap 
separates the QP peak from the rest of the spectrum at higher energy. 

Let me now study more quantitatively the behavior of the weight of the QP peak. Simi- 
larily to the one-hole case, I can define the Z-factor as the square of the overlap 
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(3.3) 



The denominator is a simple normalization factor |22j so that Zih represents the relative 



spectral weight (between and 1) located in the QP peak. 

One of the key issues of the problem of strongly correlated fermions is the possibility 



that spin and charge decoupling which takes place in one dimension [23] would also occur 



in two dimensions as proposed by Anderson |24]] and others P5| . If such a scenario is valid 
then, as a simple consequence of the spin-charge separation, one would expect the absence of 
sharp QP. It is then crutial to test this hypothesis. Previous calculations [|2]|| suggested that 
the single hole Green function in the AF background does exhibit an undumped QP pole. 
However, since this approach neglected any possible hole-hole interaction it is necessary to 
re-examine the same problem when more than a single hole is present. The case of two holes 
considered here is the first step towards the finite hole density. 

Z 2 h is plotted in Fig. lb as a function of J. Fig. lb shows that Z 2 u varies almost linearly 
with J in the range 0.1 < J < 1. It is also interesting to compare these data with the single 
hole QP weight Z\h obtained recently on the same clusters ||, 

Zih = E I 1 4,. I O I 2 , (3-4) 



y-ii - \ z 



K stands for the finite momentum of the one hole GS [17]. Table 4 shows the data for Z 2 h 
for the various clusters as well as (Z^) 2 on 26 sites. We note that, to a good approximation, 

lh) 

Although Z 2 h is quite small in the physically interesting range 0.1 < J < 0.5 only a finite 
size analysis can tell whether or not it vanishes in the bulk. The data for Z 2 h are plotted in 
Fig. 2b as a function of 1/N. It should be noted that, contrary to the one hole case P,|i~7f 
the GS momentum is identical for all the clusters and equal to 0. No simple scaling law can 
be deduced from these plots. However, the numbers at N=26 seem to indicate a convergence 
towards a finite number in the bulk. 

Lastly, let me consider the behavior of the spectral function vs the coupling constant J as 
shown in Figs. 5a-f. There are many similarities with the behavior of the single hole spectral 



function calculated on the 4x4 and larger clusters [|J. In the limit J=0 the spectrum 



in Fig. 5a is symmetric around u = 0. This is in fact easy to understand. The change 
uj — > — ijj is equivalent to t — > — t ie c ia — > —c-^ a for, let say, the even sites. Under such a 
transformation one can easily check that A* — > —A' so that eventually A 2 h{^>) = A 2 h(— u) 



for J=0. The spectral density is almost constant over an energy range of ~ 13 1. A sharp 
peak (actually a 5-function) can be seen at uj = 0. When J is turned on (Figs. 5b, c) the 
spectrum is no longer symmetric with respect to uj and no longer (rigourously) bounded 
from above (strictly speaking in the bulk) since there now exist states with excitation up to 
a magnetic energy oc N J, N being the system size. However, in pratice, the matrix elements 
in ( |3.1|) decrease exponantially fast with increasing energy outside an energy range of a few 
t (~ 10 - 12 t). As seen above a QP peak of weight oc J appears at the bottom of the band. 
However for small J most of the spectral weight lies above. The sharp uj = peak of the 
J — > limit is rapidly broaden by J. Reminiscence of this peak can still be clearly seen up to 
J=0.2. If J is increased further (figs. 5d-f) most of the features of the J = limit vanish. On 
the other hand, sharp peaks (nevertheless of finite width) appear right above the QP peak. 
These structures are characteristic of the J ~ t region and might be related to the internal 
structure of the pair. Indeed, the 2 holes (for sufficiently large J) are confined on neighboring 
sites by some effective potential which naively increases linearly with the hole separation (if 
spin fluctuations are not considered). It is interesting to note that the appearance of these 
resonances coincide with the approximate J value at which binding sets up (see Sec. II). 



IV. CONCLUSION 

In this paper, calculations of the binding energy as well as of the pair field spectral 
function (of d-wave symmetry) have been reported for the t-J model. By diagonalising 
clusters with up to 26 sites and 2 holes it has been definitely established that binding 
between holes occurs only above a critical threshold of J in agreement with previous work. 
I estimate this critical value to be in the range 0.3 - 0.5. The spectral function of the pair 
operator shows unambiguously a QP peak at the bottom of the band which does not seem 
to broaden or disappear when the system size is increased. Its weight grows almost like J 
for J < 1. In the physical region 0.1 < J < 0.5 most of the spectral weight is spread above 
the peak on a wide energy range of the order of a few t. Two qualitatively different regimes 
are found: (i) for J < 0.25, the spectral function exhibits a pronounced bump around the 
band center which eventually merges into a sharp uj = peak when J=0; (ii) for J > 0.25, 
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especially J ~ 1, sharp (small) resonances are seen above the QP peak. These structures 
might be related to the internal nature of the pair bound by some effective potential. 

The numerical calculations were done on the CRAY-2 of Centre de Calcul Vectoriel pour 
la Recherche (CCVR), Palaiseau, France. I am particularly greatful to the CCVR staff for 
helpful assistance. I also wish to thank the Supercomputer Computation Research Institute 
(Tallahassee) for hospitality at an early stage of this work and acknowledge my co-workers 
E. Dagotto, J. Riera, A. Moreo and M. Novotny for many fruitful discussions. 1 am also 
indebted to J. Riera for providing me with many data to check the computer codes. 
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TABLE IB 
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TABLE 2 
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TABLE 3 
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TABLE 4 
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TABLE CAPTIONS 



Table 1 

Symmetry properties of the half filled (A) and two hole (B) GS. The fermion convention is 
used. 

Table 2 

Two hole ground state energies of the t-J model Eq~ 2 - Eg on clusters 4x4, VT8 x Vl8, 
V20 x v 7 ^ and V^E x V^E for several values of J (0.1 to 1). 

Table 3 

Binding energies in the t-J model on clusters 4x4, VT8 x y/l8, y^20 x and \^26 x \^26 
for several values of J (0.1 to 1). Note that for 18 sites a level crossing appears in the one hole 
sector (indicated by *). 

Table 4 

Spectral weight Z 2 h of the QP peak in the pair spectral function of the t-J model on clusters 
4x4, x \/l8, V20 x \^20 and \^26 x \^26 for several values of J (0.1 to 1). For comparison, 
Zf h corresponding to the one hole GS on 26 sites is also listed on the last line (marked by *). 
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FIGURE CAPTIONS 



Figure 1 

Binding energy A 2 h (a) and spectral weight Z 2 h (b) vs J. The symbols corresponding to the 
different clusters are shown on the plot. 

Figure 2 

Binding energy (a) and spectral weight Z^h (b) plotted (arbitrarily) vs inverse system size 
1/N. The value of J corresponding to each set of points is indicated on the plot on the left side 
of the corresponding data. 

Figure 3 

Spectral function of the d-wave pair operator vs frequency (in unit of t) for N = 26 and 
J = 0.3. The Lanczos iteration number N it used to calculate the continued fraction of the 
spectral function is shown on the plots. 

Figure 4 

Spectral function of the d-wave pair operator vs frequency for increasing cluster sizes. J = 0.1 
(a - d), 0.3 (e - h) and 1 (i - 1). The frequency uj is measured in unit of t. 

Figure 5 

Spectral function of the d-wave pair operator vs frequency on the \/26 x a/26 cluster for various 
J values, uj is measured in unit of t. The continued fraction expension was truncated after 
N it = 450 iterations (except for J = where N it = 550). 
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